#ifndef GRDECLDATA_H
#define GRDECLDATA_H

#include <cstdio>
#include <string>

#include "base.h"
#include "geo2d.h"
#include "geo3d.h"

// Point usage declare here
using geo2d::Point2d;

using geo3d::Point3d;
using geo3d::Line;

struct Corn8 {
    double zval[2][2][2];
    Point3d coord[2][2][2];
    
    void print();
    Point3d GetCenter();
};

struct GridProperties {
    bool inited = false;
    bool isPillarVertical; // is all pillars vertical
    bool isGridDense; // bottom of upper grid is the same as top of lower grid
    bool isCoordMonotonic; // 是否所有格点都满足三个维度上的单调性
    int x_mono_dir, y_mono_dir, z_mono_dir; // 三个维度的单调方向

    void print();
};

struct GridRange {
    Point3d rmin, rmax;

    void print();
};

class GRDECLData {
public:
    INDEX nx, ny, nz;
    double origin_x, origin_y;
    Line* pillars;
    Corn8* corns;
    double* young;
    double* poisson;
    double* rockden;
    double* poro;
    bool* actNum; // active cells
    bool* isDense; // dense cells
    bool* isFault; // fault points
    bool* isConvex; // convex quadrilaterals
    std::vector<geo3d::Fault> faults;
    std::string filename;
    // 输出路径
    std::string output_path;

    GridProperties properties;
    GridRange range;

    double avg_young, avg_poisson, avg_rockden, avg_poro;

    GRDECLData() {}
    GRDECLData(std::string filename);

    Line& GetPillar(INDEX y, INDEX x);
    Corn8& GetCorner(INDEX z, INDEX y, INDEX x);

    double GetYoung(INDEX z, INDEX y, INDEX x);
    double GetPoisson(INDEX z, INDEX y, INDEX x);
    double GetRockDen(INDEX z, INDEX y, INDEX x);
    double GetPoro(INDEX z, INDEX y, INDEX x);

    inline Point2d GetPillarAtZval(INDEX y, INDEX x, double zval);
    Point2d GetPillarAtZval(const Line& pillar, double zval);
    inline Point2d GetPointFromPillar(INDEX y, INDEX x, double zval);

    int InitFromFile(); // 从文件中读取网格数据
    int InitFromFileOpt(); // InitFromFile 的快速读入版本
    void InitProperties(); // 初始化网格属性
    void InitRange();
    void InitConvex(); // 判断柱线cell是否为凸四边形
    void SetOutputPath(std::string output_path) { // 设置输出路径
        this->output_path = output_path;
    }

    void ReadPropsFromFileOpt(std::string filename); // ReadPropsFromFile 的快速读入版本
    void CalAverageProps(); // 计算平均属性值
    void PrintHeadProps();
    int PrintPropToFile(std::string filename);

    void ReadFaultFromFile(std::string filename);

    void WriteVTK(std::string filename, bool cellData = false);
    void WriteLayerVTK(std::string base_filename, bool cellData = false);
    void WriteFaultVTK(std::string filename);

    // 根据点坐标插值得到新柱线; new_pillar 返回新柱线; px,py 返回相邻柱线坐标
    void InterpolateNewPillar(double x, double y, double z, Line& new_pillar, INDEX& px, INDEX& py); 
    // 根据相邻柱线 px, py 计算柱线 pillar 在第 layer 层的 z 坐标
    double InterpolatePillarZvalAtLayer(int layer, Line pillar, INDEX px, INDEX py);

    // 获取点所在的单元index
    void GetCornerOfPoint(INDEX &z, INDEX &y, INDEX &x, const Point3d &p, INDEX prev_y, INDEX prev_x, std::mt19937 &gen);
    
    void GetNaiveBottomSurface(std::vector<Point3d> &pts, std::vector<INDEX> &bdrs, std::vector<INDEX> &tris);

    void GetFaultsFromFile(std::vector<Point3d> &fault_points, std::vector<INDEX> &fault_lines, std::vector<INDEX> &fault_quads);

private:
    INDEX nxy, nxyz, nx1;

    // InterpolatePillar 函数每次从 (_curr_pillar_x, _curr_pillar_y) 开始继续查找，从而利用局部性
    INDEX _curr_pillar_x = 0;
    INDEX _curr_pillar_y = 0;
};

#endif // GRDECLDATA_H
